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An self-contained elastic theory is derived which accounts both for mechanical yield and shear- 
induced volume dilatancy. Its two essential ingredients are thermodynamic instability and the 
dependence of the elastic moduli on compression. 
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Stable sand piles show clear elastic behavior, though 
conventional elasticity theory cannot be appropriate, as 
sand piles also possess a steepest slope. The associated 
angle, the Coulomb angle p, has a typical value of around 
30° for dry sand. So granular materials may be taken to 
interpolate between fluid (no elasticity, ip = 0°) and solid 
(conventional elasticity, p > 90°). This is usually under- 
stood in terms of the solid friction law: Subject to grav- 
itational pull, a sand grain resting on a slope of angle 
experiences the force N = pG cos 9 normal to the surface, 
and S = pGsinO along it, with S/N — tan 6. (p is the 
mass density and G the gravitational constant.) Since 
static friction will only sustain a maximal S/N , there is 
a maximal 0, which one may identify with the Coulomb 
angle p, or S/N < tan</?. Now, requiring this to hold 
both along any plane and everywhere in the bulk, one 
may reinterpret S and N as components of the stress ten- 
sor (Tij, respectively tangential and normal to the given 
plane. Then the inequality may, with <j\ and 03 as the 
largest and smallest Eigenvalues of <7y, be written as 

|((Ti - cr 3 )/((Ti +cr 3 )| < sva.pi. (1) 

This is the "Coulomb yield condition/ 1 or "Coulomb law 
of internal friction/ 1 a textbook formula of soil mechan- 
ics 0, employed to impose mechanical yield upon con- 
ventional elasticity. Although this formula captures es- 
sential granular physics, it possesses a number of defi- 
ciencies. First of all, the Coulomb condition "pre-empts" 
what should have been the result of a proper theory: Ide- 
ally, one would like to start from a continuum theory, 
calculate the spatial dependence of the stress with ap- 
propriate boundary conditions, and arrive - with some 
substantiated understanding - at the fact that there is 
a maximal shear stress in sand, above which the system 
is mechanically unstable. This includes especially an ex- 
pression for the Coulomb angle p. The Coulomb law 
postulates the last bit and employs it backward. 

Second, the Coulomb condition takes for granted that 
mechanical yield is determined by a unique p, indepen- 
dent of geometry. Third, it is not obvious that gran- 
ular materials behave as conventional solids up to the 
yield point, without any "precursor" behavior. Fourth, 
the Coulomb condition partially contradicts conventional 



elasticity, and our understanding is rendered regrettably 
precarious. Last not least, Reynolds dilatancy [2j - the 
volume expansion concurring with shear motion in gran- 
ular materials - should be an integral part of mechani- 
cal yield, yet is completely ignored by the Coulomb law. 
Imagine a pile of stacked steel balls, and envisage how 
a shearing displacement lift the balls from their close- 
packed positions and give rise to volume expansion 
which eventually leads to yield. In experiments Q and 
simulations this is what has been observed. 

Although the Coulomb condition appears unique to 
granular systems, its sole purpose is to account for me- 
chanical yield, or the lack of elastic configurations for 
certain values of stress and strain. Yield is a widespread 
phenomenon in many solids at high stresses, which are 
well accounted for by linear elasticity at lower values 
of stress. "Low, high" are of course relative concepts, 
and the noteworthy point is: Being characterized by a 
quadratic elastic free energy and a linear stress, 

U = \K h u\ n + K a u° ke u° ke , (2) 
<x y - = -K b u m 8ij - 2K a u\ , (3) 

linear elastic theory cannot possibly account for yield, 
as it provides stable elastic solutions for arbitrary strains 
and stresses, however high. (u nn is the trace of the strain 
tensor tiy, u kl = Uki — u nn 6ki/3 is its traceless part, and 
Kb,K a are the constant compressional and shear mod- 
uli.) Yet contrary to prevalent perception, this is not a 
general feature of elasticity: Adding nonlinear terms to 
Eqs ilL'l.'{|) may well render given elastic solutions insta- 
ble for some variable range. And it appears obvious that 
elastic instabilities, or more generally, the lack of elastic 
solutions, are to be identified with yield, the lack of elas- 
tic configurations. Doing so embraces yield as the generic 
phenomenon that it is, and does away with extraneous 
inputs such as the Coulomb condition. Note also that 
yield therefore marks the end of the range of validity for 
elasticity. Only a more comprehensive theory including 
dissipative terms is able to describe what then happens, 
typically plastic flows. 

In granular materials, one need not look far for terms 
that fit the above description. Consider two solid spheres 
in contact 5] to find U ~ /1 5 / 2 , / ~ ft 3 / 2 , where U is the 
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clastic energy, / the applied force, and h the relative 
change in height. The latter relation is not linear, be- 
cause the area of contact between the two spheres and 
the amount of compressed mass increase with h. Assum- 
ing one can scale up this two-body result to granular ma- 
terials - in usual parlance, that they possess "Hertz con- 
tacts, " we identify U, f, h respectively with f e i, a nn , u nn 
of Eqs (|2I3J) . and conclude Kb ~ ii nn 1 '' 2 ~ oVm 1 / 3 is no 
longer constant. Realistically, with more than two grains 
in contact, Kb will still depend on o~ nn and vanish with 
it, so one may more generally take Kb ~ <Xnn • And since 
the physics of increasing contact area is similar for shear, 
also take K a ~ o nr P . Evesque and de Gennes employed 
these elastic moduli in the stress of Eq J3J to successfully 
render the pressure saturating in silos |(|. 

Aiming to generalize this "quasi-elastic theory" and 
embed it into a consistent thermodynamic framework, 
we made the following observation: With Kb, K a ~ u nn a 
and the free energy / e ; retaining its form of Eq J3J , the 
stress G{j - given by general considerations essentially as 
df e i/duij - is necessarily modified: Eq © is clearly only 
correct if K a ,Kb are constant. Our serendipitous find- 
ing reported below is, the additional terms of Oij suffice 
to account for mechanical yield and volume dilatancy, 
by rendering elastic solutions instable or untenable in a 
range of parameters appropriate for granular materials. 
(Assuming as above that it is the stress which retains 
its form, no f e i exists such that cr.y = df e i/duij holds, 
because the relevant Maxwell relations are violated.) 

Basic to our approach is the assumption that a finite 
elastic region exists. This is universally accepted in soil 
mechanics, eloquently supported by de Gennes [(|, and in 
fact corroborated by the stability of sand piles. There are 
of course well argued reservations, such as those derived 
from force arches, or from the distinction between plastic 
versus elastic contacts - introduced partly due the same 
desire to understand yield 0. However, they are usu- 
ally based on an "intergranular" or microscopic point of 
view, and the connection to the elasticity of macroscopic, 
continuous media is less than clear 0. A different argu- 
ment is the possible lack of a unique displacement field 
U{. In our opinion, elasticity is, at its core, a robust the- 
ory: In spite of crystal defects, frequently rendering Ui 
ill-defined, elasticity remains valid in solid, accounting 
for its capability to sustain shear stresses - as long as the 
defects are stationary. By merely standing, sand piles 
demonstrate the same capacity, and there is no reason 
why a carefully constructed elastic theory should not be 
able to account for all its static, macroscopic behavior as 
well. After all, whether a proposed theory does exactly 
this should be its ultimate test. 

This ends the introduction. In the following, an self- 
contained clastic theory capable of accounting for me- 
chanical yield and volume dilatancy is developed. First, 
we choose an equilibrium state of arbitrary (thermody- 
namic) temperature T and packing density p c , with no 



external forces, especially gravitation, but with the atmo- 
spheric pressure present [9(- This virtual state (in spite 
of its marginalized stability and individually compressed 
grains) is taken as our system of reference, with a van- 
ishing displacement field Ui. The associated free energy 
density is fxp/m, where f\(T) is the free energy, and 
m the mass, per grain. Turning on the gravitation and 
applying external forces will further strain the granular 
material and lead to density change. If the force is small 
and applied slowly enough, this change is elastic, 

6=1- Pc/P = -U nn - (4) 

As gravitation and normal stress cram the grains, they 
lead to finite contact areas between them, and give rise 
to finite elastic moduli. The free energy density becomes 

/ = {h/m)p + \K b ul n + K a u° ke u° ke + pGz, (5) 
Kb = K b S b , K a = K a 6 a , (6) 

with K a , K b > for 6 > 0, K b , K a = for 5 < 0. (The 
particles loose contact with one another for 6 < 0.) Lin- 
ear elasticity corresponds to a, b — 0, while "Hertz con- 
tacts" imply a,b = 1/2. However, we did not find a 
watertight general reason requiring a — b, as assumed 
in the quasi- elastic theory [jj. We believe this is a ques- 
tion of clarity versus accuracy: Experiments are better 
accounted for if b is taken slightly larger than a (see be- 
low). Yet a is the rather more important exponent, which 
alone already gives rise to mechanical yield and volume 
dilatancy. So b need not be treated with the same care, 
and one may set b — a to gain great simplifications in the 
expressions - mainly because the Poisson ratio 

v = (3K h - 2K a )/(6K b + 2K a ) (7) 

then lacks critical density dependence. Finally, Eq © 
should be taken in the spirit of an expansion that holds 
only close to p c , or for <5 <C 1. 

Given Eqs ijoKil) and employing the Eulerian notation, 
2uij = \7iUj+\7 jU— ViC/feVjC/fc, the stress is determined 
by energy and momentum conservation (cf. jlfij ) to be 
o-jj = (pp - f)6ij ~ $>ij + ^ ik u k j + ^jkUki, where *y = 
(df/duij) PtT and p = (df/dp) UtT - Keeping only the 
dominant of the nonlinear terms, we have 

aij = -K b u nn 6ij - 2K a u° j (8) 
+ 6- 1 {\bK b u 2 nn + aK a u° ek u° ik )6 tJ . 

Setting a, b = 0, only the first line, or the expression of 
linear elasticity, remains. The next line comes either from 
deriving K a , Kb with respect to p = p c (l—6), leading to a 
contribution in the chemical potential p, or equivalently, 
with respect to u nn — —6, giving rise to additional terms 
in ^ij. Although a = b ^ in the quasi-elastic theory 0], 
neither does its stress contain the second line. As will be 
shown in a future work, this is qualitatively alright for 
silos, because yield is never a problem here. 
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of sand subject to a hydrostatic pressure p and a devi- 
atory normal stress q; (iii) sand on a slope: an infinite 
layer of sand on a rough, inclined plane with the angle 
9. The implication of Eq for these experiments are, 
respectively 

S/N<twnpr, q/(2p + q) < smip 2 ; < ip 1: (11) 

with (pi — if2 denoting the Coulomb angle. Using the 
above granular elastic theory, setting for simplicity a = b, 
the respective results may also be dressed as the same 
inequalities, though (fi,if2 are now explicitly given, 



FIG. 1: Compression 8 versus pressure P (left) and shear u s 
(right). Stability is maintained only where the lines are solid. 



To obtain a feeling for the implications of Eq JHJ, con- 
sider the "pressure" P = <7fcfc/3, as a function of the 
compression 5 and the shear u s 



P = (1 + b/2)K b 6 1+b + akaul/5 1 - 



(9) 



8 is plotted versus P for given u s and realistic values of 
a, b, in the left half of Fig.l. The solid part of the lines 
are the stable, physical region, with a positive compressi- 
bility; the dashed lines, showing P — > 00 for 6 — > are 
the unstable, unphysical region. The reason for this can 
be understood from the right half of Fig.l, a plot of 5 
versus u s for given P. Starting at the top, with finite 
compression and no shear, 0" = 4 x 10~ 4 , u s = 0, we 
see how the solid lines decrease, or how volume dilates, 
for increasing shear. (Conventional linear elasticity or 
the elasto-plastic theory yield a straight horizontal 
line, stopping at the value asserted by the Coulomb yield 
condition.) Shear values right of the parabola-like curves 
obviously do not have elastic solutions, though as men- 
tioned plastic flow solutions of course must exist. 

Thermodynamic stability, however, is lost even before 
the vanishing of the elastic solutions, at the points given 
by the arrows, where the solid lines turn into dashed ones. 
This is because / of Eq (JSJ is convex only for certain 
values of S and u s : With / = \K b 5 h+2 + K a S a u 2 + • • • , 
thermodynamic stability requires (d 2 f /dS 2 )(d 2 f /du 2 ) > 
(d 2 f/dSdu s ) 2 , or 



u 2 JS 2 < (2 + 6)(1 + b)K b /[2a(l + a)K a ] 



Dashed lines of both figures represent parameters that 
do not satisfy this condition. Note that if the medium is 
not at all compressed, no finite shear is stable. 

For further comparison of this theory to the Coulomb 
yield condition (JTJ, we shall in the following consider 
three typical experimental setups (see upper inset of 
Fig. 2): (i) simple shear test: an infinite layer of sand 
subject to a normal and a shear force density, N and 
S; (ii) axisymmetric triaxial test: a cylindrical sample 
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These expressions are easily derived. For the simple shear 
test, the symmetry of the geometry and the force balance 
VfeCTifc = (neglecting gravity) prescribe constant strain, 
with the displacement given as U x , U y ~ y, U z = 0. In- 
serting the nonvanishing components of the strain, u xy 
and u yy = — S, into Eq JSJ, we obtain 



= (1 + \)K t S + |(1 + %)K a S + 2aK a u 2 /5, (14) 



(J X y = ~2K a u 



■>"y 



®~ XX = O-yy - 2K a S. (15) 



The boundary conditions impose o~ y 



N, 



S. 



Solving 8 and u xy as functions of N and S, we find that 
no solution exists if N/S exceeds a maximal value. In 
addition, S and u xy must satisfy Eq Hl()[) . which for the 
present case is the more stringent one, leading to Eq l|12|l . 
Note since of Eq H14[) is similar to Eq © , (they are struc- 
turally identical for a — 6), it also displays dilatancy. 

The strain is again constant for the axisymmetric tri- 
axial test. With the displacement vector given as U x ~ x, 
U y ~ y, U z ~ z, the strain is u xx — u yy = — (5 + 
AH/H )/2, u zz = AH/H Q , and u tj = for i ^ j. (H is 
the height for p = q = 0, and AH its change.) Inserting 
these into Eq JSJ) yields the stress tensor, 



j yy - %q 2 /{K a 5)-\q + K b (l 
7 ZZ = <T yy -(6 + 3AH/H a )K a . 



)S, 



(16) 



The boundary conditions are: a xx — a yy = p, o~ zz = p+q, 
(10) — f° r i 7^ j- Again, 6 and AH do not have solutions 



if q/(2p + q) exceeds the value given by Eq JEJ. [Eq flT)|) 

yields exactly the same constraint here.] 

Due to high external forces, the actual deformations 
in both above experiments tend to contain considerable 
plastic contributions, invalidating a direct comparison 
with the present theory. This is not the case for the third 
experiment, sand on a slope, as its deformation is due to 
the comparatively small gravity. The medium is uniform 
along the x and z directions (see inset of Fig. 2), so the 
displacement is of the form U x = U x {y), U y — U y (y), 
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FIG. 2: Plots of ifi (full lines) and ip2 — <pi (dashed lines) as 
functions of a for v =0.2, 0.3, 0.4 (denoted as 1, 2, 3). 



U z = 0. The stress tensor is then given by Eq JSJ and 
by integrating da yy /dy — — pG cos 9, da xy /dy = pGsmO 
under the constraint that the stress vanishes at the free 
surface (y = H). 

<J XX = (l + ±)K b 5 + l(a-l)K a 6 + 2aK a u 2 xy /6, 

a xy = -2pu xy = -GM(y) sin 9, 

<Jyy = a xx + 2K a S = GM{y) cos 9. (17) 

M(y) = j H p(y')dy' denotes mass per unit area between 
y and the free surface H; the total mass M(y — 0) is 
taken as a constant that does not vary with 9, so both 
p and H depend on 9. As in the simple shear test, the 
stability condition Eq (|10fl is the more stringent one con- 
fining the value of sin 9 and leads to Eqs i|12|) . 

Finally, we consider what these results imply for the 
parameters of the above granular elastic theory, espe- 
cially the powers a, 6 and the Poisson ratio v. In Fig. 2, 
ipi and <f2 — fi are plotted as functions of a with differ- 
ent values of v, rendered in full and dashed lines, respec- 
tively. The angle ip% varies between 15° and 35°, with 
1^2 — <pi \ /f i 5% in the parameter regime 0.2 < a < 1 
and 0.35 < v < 0.4. Insisting on tf2 = <fi = 30° for 
Coulomb media, we find a — 0.27 and v — 0.36. (The 
fact that |(^2 — diverges for small a indicates that lin- 
ear elasticity, a — > 0, is incompatible with a unique angle 
of friction.) Because (pi is the maximum angle for in- 
clined planes with a layer of sand at rest, we may identify 
(pi with Bangold's maximum angle of stability. Experi- 
mentally, this angle is found to vary from 30° (spherical 
grains) to 60 — 70° (angular and rough grains) [12J . Ac- 
cording to Fig. 2, this implies an a between and 0.4. 

Because we assumed b = a, b also varies between 
and 0.4. With P ~ S 1+b or K b ~ P 13 , (3 = 6/(1 + b), cf. 
Eqs © and JJjJ, we have (5 varying between and 0.28 
- smaller than the i of the Hertz contact 0. This is an 
indication that a — b oversimplifies. Repeating the above 



calculation allowing a ^ b, we find <p\ = <p% = 30° by 
taking a = 0.4, b = 0.5 (ie. (3 = §), and K a /K b = 0.36. 

At higher pressures (from 700 to 7000 kg/m 2 ), = 0.5 
(or 6=1) was measured and referred to as the P 1 / 2 - 
dependence |13|). Various microscopic reasons were pro- 
posed for this deviation from the Hertz law 0, Il5| . 
Within the present framework, this is easily accounted 
for by the packing dependence of b(p c ), a feature that we 
shall study in future works. 

It is important to realize that the Poisson ratio v as 
given in Eq J7J is, in granular materials, not the same as 
that given by 

_ ( du xx ldg \ = 1 6 2 + 36+(6 2 + 36 + 6)t/ 
Vtri \du zz /dq) 2 6 2 + 36 + 3 + 6(6 + 3)^ 1 ' 

(at q — > 0), eg measured in triaxial experiments. [The 
second expression is obtained by employing Eq (|16|l .] 
Taking 6 = 0.5, a — 0.4, K a /K b = 0.36 as above, and 
S = 10~ 5 , we obtain v tr i — 0.25, well comparable to the 
measured value of around 0.17 to 0.25 [lfj. (Enforcing 
a = b again leads to discrepancy: Taking a — b — 0.27 
and ^=0.36 as given above, we have vtri = 0.4.) 
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